!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!

! **************************************************************************************************
!> \brief A collection of methods to treat the QM/MM electrostatic coupling
!> \par History
!>      5.2004 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
MODULE qmmm_gpw_energy
   USE cell_types,                      ONLY: cell_type,&
                                              pbc
   USE cp_control_types,                ONLY: dft_control_type
   USE cp_log_handling,                 ONLY: cp_get_default_logger,&
                                              cp_logger_type
   USE cp_output_handling,              ONLY: cp_p_file,&
                                              cp_print_key_finished_output,&
                                              cp_print_key_should_output,&
                                              cp_print_key_unit_nr
   USE cp_realspace_grid_cube,          ONLY: cp_pw_to_cube
   USE cp_spline_utils,                 ONLY: pw_prolongate_s3,&
                                              spline3_nopbc_interp,&
                                              spline3_pbc_interp
   USE cube_utils,                      ONLY: cube_info_type
   USE input_constants,                 ONLY: do_par_atom,&
                                              do_qmmm_coulomb,&
                                              do_qmmm_gauss,&
                                              do_qmmm_none,&
                                              do_qmmm_pcharge,&
                                              do_qmmm_swave
   USE input_section_types,             ONLY: section_get_ivals,&
                                              section_vals_get_subs_vals,&
                                              section_vals_type,&
                                              section_vals_val_get
   USE kinds,                           ONLY: dp
   USE message_passing,                 ONLY: mp_para_env_type
   USE mm_collocate_potential,          ONLY: collocate_gf_rspace_NoPBC
   USE particle_list_types,             ONLY: particle_list_type
   USE particle_types,                  ONLY: particle_type
   USE pw_env_types,                    ONLY: pw_env_get,&
                                              pw_env_type
   USE pw_methods,                      ONLY: pw_zero
   USE pw_pool_types,                   ONLY: pw_pool_p_type,&
                                              pw_pools_create_pws,&
                                              pw_pools_give_back_pws
   USE pw_types,                        ONLY: pw_r3d_rs_type
   USE qmmm_gaussian_types,             ONLY: qmmm_gaussian_p_type,&
                                              qmmm_gaussian_type
   USE qmmm_se_energy,                  ONLY: build_se_qmmm_matrix
   USE qmmm_tb_methods,                 ONLY: build_tb_qmmm_matrix,&
                                              build_tb_qmmm_matrix_pc,&
                                              build_tb_qmmm_matrix_zero
   USE qmmm_types_low,                  ONLY: qmmm_env_qm_type,&
                                              qmmm_per_pot_p_type,&
                                              qmmm_per_pot_type,&
                                              qmmm_pot_p_type,&
                                              qmmm_pot_type
   USE qmmm_util,                       ONLY: spherical_cutoff_factor
   USE qs_environment_types,            ONLY: get_qs_env,&
                                              qs_environment_type
   USE qs_ks_qmmm_types,                ONLY: qs_ks_qmmm_env_type
   USE qs_subsys_types,                 ONLY: qs_subsys_get,&
                                              qs_subsys_type
#include "./base/base_uses.f90"

   IMPLICIT NONE
   PRIVATE

   LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_gpw_energy'

   PUBLIC :: qmmm_el_coupling
   PUBLIC :: qmmm_elec_with_gaussian, &
             qmmm_elec_with_gaussian_LR, qmmm_elec_with_gaussian_LG
!***
CONTAINS

! **************************************************************************************************
!> \brief Main Driver to compute the QM/MM Electrostatic Coupling
!> \param qs_env ...
!> \param qmmm_env ...
!> \param mm_particles ...
!> \param mm_cell ...
!> \par History
!>      05.2004 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   SUBROUTINE qmmm_el_coupling(qs_env, qmmm_env, mm_particles, mm_cell)
      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(qmmm_env_qm_type), POINTER                    :: qmmm_env
      TYPE(particle_type), DIMENSION(:), POINTER         :: mm_particles
      TYPE(cell_type), POINTER                           :: mm_cell

      CHARACTER(len=*), PARAMETER                        :: routineN = 'qmmm_el_coupling'

      INTEGER                                            :: handle, iw, iw2
      LOGICAL                                            :: mpi_io
      TYPE(cp_logger_type), POINTER                      :: logger
      TYPE(dft_control_type), POINTER                    :: dft_control
      TYPE(mp_para_env_type), POINTER                    :: para_env
      TYPE(particle_list_type), POINTER                  :: particles
      TYPE(pw_env_type), POINTER                         :: pw_env
      TYPE(pw_pool_p_type), DIMENSION(:), POINTER        :: pw_pools
      TYPE(qs_ks_qmmm_env_type), POINTER                 :: ks_qmmm_env_loc
      TYPE(qs_subsys_type), POINTER                      :: subsys
      TYPE(section_vals_type), POINTER                   :: input_section, interp_section, &
                                                            print_section

      CALL timeset(routineN, handle)
      logger => cp_get_default_logger()
      NULLIFY (ks_qmmm_env_loc, pw_pools, pw_env, input_section, dft_control)
      CALL get_qs_env(qs_env=qs_env, &
                      pw_env=pw_env, &
                      para_env=para_env, &
                      input=input_section, &
                      ks_qmmm_env=ks_qmmm_env_loc, &
                      subsys=subsys, &
                      dft_control=dft_control)
      CALL qs_subsys_get(subsys, particles=particles)

      CALL pw_env_get(pw_env=pw_env, pw_pools=pw_pools)
      print_section => section_vals_get_subs_vals(input_section, "QMMM%PRINT")
      iw = cp_print_key_unit_nr(logger, print_section, "PROGRAM_RUN_INFO", &
                                extension=".qmmmLog")
      IF (iw > 0) &
         WRITE (iw, '(T2,"QMMM|",1X,A)') "Information on the QM/MM Electrostatic Potential:"
      !
      ! Initializing vectors:
      !        Zeroing v_qmmm_rspace
      CALL pw_zero(ks_qmmm_env_loc%v_qmmm_rspace)
      IF (dft_control%qs_control%semi_empirical) THEN
         ! SEMIEMPIRICAL
         SELECT CASE (qmmm_env%qmmm_coupl_type)
         CASE (do_qmmm_coulomb, do_qmmm_none)
            CALL build_se_qmmm_matrix(qs_env, qmmm_env, mm_particles, mm_cell, para_env)
            IF (qmmm_env%qmmm_coupl_type == do_qmmm_none) THEN
               IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') &
                  "No QM/MM Electrostatic coupling. Just Mechanical Coupling!"
            END IF
         CASE (do_qmmm_pcharge)
            CPABORT("Point charge QM/MM electrostatic coupling not yet implemented for SE.")
         CASE (do_qmmm_gauss, do_qmmm_swave)
            CPABORT("GAUSS or SWAVE QM/MM electrostatic coupling not yet implemented for SE.")
         CASE DEFAULT
            IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') "Unknown Coupling..."
            CPABORT("")
         END SELECT
      ELSEIF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
         ! DFTB
         SELECT CASE (qmmm_env%qmmm_coupl_type)
         CASE (do_qmmm_none)
            IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') &
               "No QM/MM Electrostatic coupling. Just Mechanical Coupling!"
            CALL build_tb_qmmm_matrix_zero(qs_env, para_env)
         CASE (do_qmmm_coulomb)
            CALL build_tb_qmmm_matrix(qs_env, qmmm_env, mm_particles, mm_cell, para_env)
         CASE (do_qmmm_pcharge)
            CALL build_tb_qmmm_matrix_pc(qs_env, qmmm_env, mm_particles, mm_cell, para_env)
         CASE (do_qmmm_gauss, do_qmmm_swave)
            CPABORT("GAUSS or SWAVE QM/MM electrostatic coupling not implemented for DFTB.")
         CASE DEFAULT
            IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') "Unknown Coupling..."
            CPABORT("")
         END SELECT
      ELSE
         ! QS
         SELECT CASE (qmmm_env%qmmm_coupl_type)
         CASE (do_qmmm_coulomb)
            CPABORT("Coulomb QM/MM electrostatic coupling not implemented for GPW/GAPW.")
         CASE (do_qmmm_pcharge)
            CPABORT("Point Charge QM/MM electrostatic coupling not implemented for GPW/GAPW.")
         CASE (do_qmmm_gauss, do_qmmm_swave)
            IF (iw > 0) &
               WRITE (iw, '(T2,"QMMM|",1X,A)') &
               "QM/MM Coupling computed collocating the Gaussian Potential Functions."
            interp_section => section_vals_get_subs_vals(input_section, &
                                                         "QMMM%INTERPOLATOR")
            CALL qmmm_elec_with_gaussian(qmmm_env=qmmm_env, &
                                         v_qmmm=ks_qmmm_env_loc%v_qmmm_rspace, &
                                         mm_particles=mm_particles, &
                                         aug_pools=qmmm_env%aug_pools, &
                                         para_env=para_env, &
                                         eps_mm_rspace=qmmm_env%eps_mm_rspace, &
                                         cube_info=ks_qmmm_env_loc%cube_info, &
                                         pw_pools=pw_pools, &
                                         auxbas_grid=qmmm_env%gridlevel_info%auxbas_grid, &
                                         coarser_grid=qmmm_env%gridlevel_info%coarser_grid, &
                                         interp_section=interp_section, &
                                         mm_cell=mm_cell)
         CASE (do_qmmm_none)
            IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') &
               "No QM/MM Electrostatic coupling. Just Mechanical Coupling!"
         CASE DEFAULT
            IF (iw > 0) WRITE (iw, '(T2,"QMMM|",1X,A)') "Unknown Coupling..."
            CPABORT("")
         END SELECT
         ! Dump info on the electrostatic potential if requested
         IF (BTEST(cp_print_key_should_output(logger%iter_info, print_section, &
                                              "POTENTIAL"), cp_p_file)) THEN
            mpi_io = .TRUE.
            iw2 = cp_print_key_unit_nr(logger, print_section, "POTENTIAL", &
                                       extension=".qmmmLog", mpi_io=mpi_io)
            CALL cp_pw_to_cube(ks_qmmm_env_loc%v_qmmm_rspace, iw2, &
                               particles=particles, &
                               stride=section_get_ivals(print_section, "POTENTIAL%STRIDE"), &
                               title="QM/MM: MM ELECTROSTATIC POTENTIAL ", &
                               mpi_io=mpi_io)
            CALL cp_print_key_finished_output(iw2, logger, print_section, &
                                              "POTENTIAL", mpi_io=mpi_io)
         END IF
      END IF
      CALL cp_print_key_finished_output(iw, logger, print_section, &
                                        "PROGRAM_RUN_INFO")
      CALL timestop(handle)
   END SUBROUTINE qmmm_el_coupling

! **************************************************************************************************
!> \brief Compute the QM/MM electrostatic Interaction collocating the gaussian
!>      Electrostatic Potential
!> \param qmmm_env ...
!> \param v_qmmm ...
!> \param mm_particles ...
!> \param aug_pools ...
!> \param cube_info ...
!> \param para_env ...
!> \param eps_mm_rspace ...
!> \param pw_pools ...
!> \param auxbas_grid ...
!> \param coarser_grid ...
!> \param interp_section ...
!> \param mm_cell ...
!> \par History
!>      06.2004 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   SUBROUTINE qmmm_elec_with_gaussian(qmmm_env, v_qmmm, mm_particles, &
                                      aug_pools, cube_info, para_env, eps_mm_rspace, pw_pools, &
                                      auxbas_grid, coarser_grid, interp_section, mm_cell)
      TYPE(qmmm_env_qm_type), POINTER                    :: qmmm_env
      TYPE(pw_r3d_rs_type), INTENT(IN)                   :: v_qmmm
      TYPE(particle_type), DIMENSION(:), POINTER         :: mm_particles
      TYPE(pw_pool_p_type), DIMENSION(:), POINTER        :: aug_pools
      TYPE(cube_info_type), DIMENSION(:), POINTER        :: cube_info
      TYPE(mp_para_env_type), POINTER                    :: para_env
      REAL(KIND=dp), INTENT(IN)                          :: eps_mm_rspace
      TYPE(pw_pool_p_type), DIMENSION(:), POINTER        :: pw_pools
      INTEGER, INTENT(IN)                                :: auxbas_grid, coarser_grid
      TYPE(section_vals_type), POINTER                   :: interp_section
      TYPE(cell_type), POINTER                           :: mm_cell

      CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_elec_with_gaussian'

      INTEGER                                            :: handle, handle2, igrid, ilevel, &
                                                            kind_interp, lb(3), ngrids, ub(3)
      LOGICAL                                            :: shells
      TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:)    :: grids

      CPASSERT(ASSOCIATED(mm_particles))
      CPASSERT(ASSOCIATED(qmmm_env%mm_atom_chrg))
      CPASSERT(ASSOCIATED(qmmm_env%mm_atom_index))
      CPASSERT(ASSOCIATED(aug_pools))
      CPASSERT(ASSOCIATED(pw_pools))
      !Statements
      CALL timeset(routineN, handle)
      ngrids = SIZE(pw_pools)
      CALL pw_pools_create_pws(aug_pools, grids)
      DO igrid = 1, ngrids
         CALL pw_zero(grids(igrid))
      END DO

      shells = .FALSE.

      CALL qmmm_elec_with_gaussian_low(grids, mm_particles, &
                                       qmmm_env%mm_atom_chrg, qmmm_env%mm_atom_index, &
                                       cube_info, para_env, eps_mm_rspace, qmmm_env%pgfs, &
                                       auxbas_grid, coarser_grid, qmmm_env%potentials, &
                                       mm_cell=mm_cell, dOmmOqm=qmmm_env%dOmmOqm, periodic=qmmm_env%periodic, &
                                       per_potentials=qmmm_env%per_potentials, par_scheme=qmmm_env%par_scheme, &
                                       qmmm_spherical_cutoff=qmmm_env%spherical_cutoff, shells=shells)

      IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
         CALL qmmm_elec_with_gaussian_low(grids, qmmm_env%added_charges%added_particles, &
                                          qmmm_env%added_charges%mm_atom_chrg, &
                                          qmmm_env%added_charges%mm_atom_index, &
                                          cube_info, para_env, eps_mm_rspace, qmmm_env%added_charges%pgfs, auxbas_grid, &
                                          coarser_grid, qmmm_env%added_charges%potentials, &
                                          mm_cell=mm_cell, dOmmOqm=qmmm_env%dOmmOqm, periodic=qmmm_env%periodic, &
                                          per_potentials=qmmm_env%added_charges%per_potentials, par_scheme=qmmm_env%par_scheme, &
                                          qmmm_spherical_cutoff=qmmm_env%spherical_cutoff, shells=shells)
      END IF
      IF (qmmm_env%added_shells%num_mm_atoms > 0) THEN
         shells = .TRUE.
         CALL qmmm_elec_with_gaussian_low(grids, qmmm_env%added_shells%added_particles, &
                                          qmmm_env%added_shells%mm_core_chrg, &
                                          qmmm_env%added_shells%mm_core_index, &
                                          cube_info, para_env, eps_mm_rspace, qmmm_env%added_shells%pgfs, auxbas_grid, &
                                          coarser_grid, qmmm_env%added_shells%potentials, &
                                          mm_cell=mm_cell, dOmmOqm=qmmm_env%dOmmOqm, periodic=qmmm_env%periodic, &
                                          per_potentials=qmmm_env%added_shells%per_potentials, &
                                          par_scheme=qmmm_env%par_scheme, qmmm_spherical_cutoff=qmmm_env%spherical_cutoff, &
                                          shells=shells)
      END IF
      ! Sumup all contributions according the parallelization scheme
      IF (qmmm_env%par_scheme == do_par_atom) THEN
         DO ilevel = 1, SIZE(grids)
            CALL para_env%sum(grids(ilevel)%array)
         END DO
      END IF
      ! RealSpace Interpolation
      CALL section_vals_val_get(interp_section, "kind", i_val=kind_interp)
      SELECT CASE (kind_interp)
      CASE (spline3_nopbc_interp, spline3_pbc_interp)
         ! Spline Iterpolator
         CALL para_env%sync()
         CALL timeset(TRIM(routineN)//":spline3Int", handle2)
         DO Ilevel = coarser_grid, auxbas_grid + 1, -1
            CALL pw_prolongate_s3(grids(Ilevel), &
                                  grids(Ilevel - 1), &
                                  aug_pools(Ilevel)%pool, &
                                  param_section=interp_section)
         END DO
         CALL timestop(handle2)
      CASE DEFAULT
         CPABORT("")
      END SELECT
      lb = v_qmmm%pw_grid%bounds_local(1, :)
      ub = v_qmmm%pw_grid%bounds_local(2, :)

      v_qmmm%array = grids(auxbas_grid)%array(lb(1):ub(1), &
                                              lb(2):ub(2), &
                                              lb(3):ub(3))

      CALL pw_pools_give_back_pws(aug_pools, grids)

      CALL timestop(handle)
   END SUBROUTINE qmmm_elec_with_gaussian

! **************************************************************************************************
!> \brief Compute the QM/MM electrostatic Interaction collocating the gaussian
!>      Electrostatic Potential - Low Level
!> \param tmp_grid ...
!> \param mm_particles ...
!> \param mm_charges ...
!> \param mm_atom_index ...
!> \param cube_info ...
!> \param para_env ...
!> \param eps_mm_rspace ...
!> \param pgfs ...
!> \param auxbas_grid ...
!> \param coarser_grid ...
!> \param potentials ...
!> \param mm_cell ...
!> \param dOmmOqm ...
!> \param periodic ...
!> \param per_potentials ...
!> \param par_scheme ...
!> \param qmmm_spherical_cutoff ...
!> \param shells ...
!> \par History
!>      06.2004 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   SUBROUTINE qmmm_elec_with_gaussian_low(tmp_grid, mm_particles, mm_charges, &
                                          mm_atom_index, cube_info, para_env, &
                                          eps_mm_rspace, pgfs, auxbas_grid, coarser_grid, &
                                          potentials, mm_cell, dOmmOqm, periodic, per_potentials, par_scheme, &
                                          qmmm_spherical_cutoff, shells)
      TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN)     :: tmp_grid
      TYPE(particle_type), DIMENSION(:), POINTER         :: mm_particles
      REAL(KIND=dp), DIMENSION(:), POINTER               :: mm_charges
      INTEGER, DIMENSION(:), POINTER                     :: mm_atom_index
      TYPE(cube_info_type), DIMENSION(:), POINTER        :: cube_info
      TYPE(mp_para_env_type), POINTER                    :: para_env
      REAL(KIND=dp), INTENT(IN)                          :: eps_mm_rspace
      TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER  :: pgfs
      INTEGER, INTENT(IN)                                :: auxbas_grid, coarser_grid
      TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER       :: potentials
      TYPE(cell_type), POINTER                           :: mm_cell
      REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: dOmmOqm
      LOGICAL, INTENT(IN)                                :: periodic
      TYPE(qmmm_per_pot_p_type), DIMENSION(:), POINTER   :: per_potentials
      INTEGER, INTENT(IN)                                :: par_scheme
      REAL(KIND=dp), INTENT(IN)                          :: qmmm_spherical_cutoff(2)
      LOGICAL, INTENT(IN)                                :: shells

      CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_elec_with_gaussian_low', &
         routineNb = 'qmmm_elec_gaussian_low'

      INTEGER                                            :: handle, handle2, IGauss, ilevel, Imm, &
                                                            IndMM, IRadTyp, LIndMM, myind, &
                                                            n_rep_real(3)
      INTEGER, DIMENSION(2, 3)                           :: bo2
      REAL(KIND=dp)                                      :: alpha, height, sph_chrg_factor, W
      REAL(KIND=dp), DIMENSION(3)                        :: ra
      REAL(KIND=dp), DIMENSION(:), POINTER               :: xdat, ydat, zdat
      TYPE(qmmm_gaussian_type), POINTER                  :: pgf
      TYPE(qmmm_per_pot_type), POINTER                   :: per_pot
      TYPE(qmmm_pot_type), POINTER                       :: pot

      NULLIFY (pgf, pot, per_pot, xdat, ydat, zdat)
      CALL timeset(routineN, handle)
      CALL timeset(routineNb//"_G", handle2)
      bo2 = tmp_grid(auxbas_grid)%pw_grid%bounds
      ALLOCATE (xdat(bo2(1, 1):bo2(2, 1)))
      ALLOCATE (ydat(bo2(1, 2):bo2(2, 2)))
      ALLOCATE (zdat(bo2(1, 3):bo2(2, 3)))
      IF (par_scheme == do_par_atom) myind = 0
      Radius: DO IRadTyp = 1, SIZE(pgfs)
         pgf => pgfs(IRadTyp)%pgf
         pot => potentials(IRadTyp)%pot
         n_rep_real = 0
         IF (periodic) THEN
            per_pot => per_potentials(IRadTyp)%pot
            n_rep_real = per_pot%n_rep_real
         END IF
         Gaussian: DO IGauss = 1, pgf%Number_of_Gaussians
            alpha = 1.0_dp/pgf%Gk(IGauss)
            alpha = alpha*alpha
            height = pgf%Ak(IGauss)
            ilevel = pgf%grid_level(IGauss)
            Atoms: DO Imm = 1, SIZE(pot%mm_atom_index)
               IF (par_scheme == do_par_atom) THEN
                  myind = myind + 1
                  IF (MOD(myind, para_env%num_pe) /= para_env%mepos) CYCLE Atoms
               END IF
               LIndMM = pot%mm_atom_index(Imm)
               IndMM = mm_atom_index(LIndMM)
               IF (shells) THEN
                  ra(:) = pbc(mm_particles(LIndMM)%r - dOmmOqm, mm_cell) + dOmmOqm
               ELSE
                  ra(:) = pbc(mm_particles(IndMM)%r - dOmmOqm, mm_cell) + dOmmOqm
               END IF
               W = mm_charges(LIndMM)*height
               ! Possible Spherical Cutoff
               IF (qmmm_spherical_cutoff(1) > 0.0_dp) THEN
                  CALL spherical_cutoff_factor(qmmm_spherical_cutoff, ra, sph_chrg_factor)
                  W = W*sph_chrg_factor
               END IF
               IF (ABS(W) <= EPSILON(0.0_dp)) CYCLE Atoms
               CALL collocate_gf_rspace_NoPBC(zetp=alpha, &
                                              rp=ra, &
                                              scale=-1.0_dp, &
                                              W=W, &
                                              pwgrid=tmp_grid(ilevel), &
                                              cube_info=cube_info(ilevel), &
                                              eps_mm_rspace=eps_mm_rspace, &
                                              xdat=xdat, &
                                              ydat=ydat, &
                                              zdat=zdat, &
                                              bo2=bo2, &
                                              n_rep_real=n_rep_real, &
                                              mm_cell=mm_cell)
            END DO Atoms
         END DO Gaussian
      END DO Radius
      IF (ASSOCIATED(xdat)) THEN
         DEALLOCATE (xdat)
      END IF
      IF (ASSOCIATED(ydat)) THEN
         DEALLOCATE (ydat)
      END IF
      IF (ASSOCIATED(zdat)) THEN
         DEALLOCATE (zdat)
      END IF
      CALL timestop(handle2)
      CALL timeset(routineNb//"_R", handle2)
      IF (periodic) THEN
         ! Long Range Part of the QM/MM Potential with Gaussians With Periodic Boundary Conditions
         CALL qmmm_elec_with_gaussian_LG(pgfs=pgfs, &
                                         cgrid=tmp_grid(coarser_grid), &
                                         mm_charges=mm_charges, &
                                         mm_atom_index=mm_atom_index, &
                                         mm_particles=mm_particles, &
                                         para_env=para_env, &
                                         per_potentials=per_potentials, &
                                         mm_cell=mm_cell, &
                                         dOmmOqm=dOmmOqm, &
                                         par_scheme=par_scheme, &
                                         qmmm_spherical_cutoff=qmmm_spherical_cutoff, &
                                         shells=shells)
      ELSE
         ! Long Range Part of the QM/MM Potential with Gaussians
         CALL qmmm_elec_with_gaussian_LR(pgfs=pgfs, &
                                         grid=tmp_grid(coarser_grid), &
                                         mm_charges=mm_charges, &
                                         mm_atom_index=mm_atom_index, &
                                         mm_particles=mm_particles, &
                                         para_env=para_env, &
                                         potentials=potentials, &
                                         mm_cell=mm_cell, &
                                         dOmmOqm=dOmmOqm, &
                                         par_scheme=par_scheme, &
                                         qmmm_spherical_cutoff=qmmm_spherical_cutoff, &
                                         shells=shells)
      END IF
      CALL timestop(handle2)
      CALL timestop(handle)

   END SUBROUTINE qmmm_elec_with_gaussian_low

! **************************************************************************************************
!> \brief Compute the QM/MM electrostatic Interaction collocating
!>      (1/R - Sum_NG Gaussians) on the coarser grid level in G-SPACE
!>      Long Range QM/MM Electrostatic Potential with Gaussian - Low Level
!>      PERIODIC BOUNDARY CONDITION VERSION
!> \param pgfs ...
!> \param cgrid ...
!> \param mm_charges ...
!> \param mm_atom_index ...
!> \param mm_particles ...
!> \param para_env ...
!> \param per_potentials ...
!> \param mm_cell ...
!> \param dOmmOqm ...
!> \param par_scheme ...
!> \param qmmm_spherical_cutoff ...
!> \param shells ...
!> \par History
!>      07.2004 created [tlaino]
!> \author Teodoro Laino
!> \note
!>      This version includes the explicit code of Eval_Interp_Spl3_pbc
!>      in order to achieve better performance
! **************************************************************************************************
   SUBROUTINE qmmm_elec_with_gaussian_LG(pgfs, cgrid, mm_charges, mm_atom_index, &
                                         mm_particles, para_env, per_potentials, &
                                         mm_cell, dOmmOqm, par_scheme, qmmm_spherical_cutoff, shells)
      TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER  :: pgfs
      TYPE(pw_r3d_rs_type), INTENT(IN)                   :: cgrid
      REAL(KIND=dp), DIMENSION(:), POINTER               :: mm_charges
      INTEGER, DIMENSION(:), POINTER                     :: mm_atom_index
      TYPE(particle_type), DIMENSION(:), POINTER         :: mm_particles
      TYPE(mp_para_env_type), POINTER                    :: para_env
      TYPE(qmmm_per_pot_p_type), DIMENSION(:), POINTER   :: per_potentials
      TYPE(cell_type), POINTER                           :: mm_cell
      REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: dOmmOqm
      INTEGER, INTENT(IN)                                :: par_scheme
      REAL(KIND=dp), DIMENSION(2), INTENT(IN)            :: qmmm_spherical_cutoff
      LOGICAL                                            :: shells

      CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_elec_with_gaussian_LG'

      INTEGER                                            :: handle, i, ii1, ii2, ii3, ii4, ij1, ij2, &
                                                            ij3, ij4, ik1, ik2, ik3, ik4, Imm, &
                                                            IndMM, IRadTyp, ivec(3), j, k, LIndMM, &
                                                            my_j, my_k, myind, npts(3)
      INTEGER, DIMENSION(2, 3)                           :: bo, gbo
      REAL(KIND=dp) :: a1, a2, a3, abc_X(4, 4), abc_X_Y(4), b1, b2, b3, c1, c2, c3, d1, d2, d3, &
         dr1, dr1c, dr2, dr2c, dr3, dr3c, e1, e2, e3, f1, f2, f3, g1, g2, g3, h1, h2, h3, p1, p2, &
         p3, q1, q2, q3, qt, r1, r2, r3, rt1, rt2, rt3, rv1, rv2, rv3, s1, s2, s3, s4, &
         sph_chrg_factor, t1, t2, t3, t4, u1, u2, u3, v1, v2, v3, v4, val, xd1, xd2, xd3, xs1, &
         xs2, xs3
      REAL(KIND=dp), DIMENSION(3)                        :: ra, vec
      REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: grid, grid2
      TYPE(pw_r3d_rs_type), POINTER                      :: pw
      TYPE(qmmm_per_pot_type), POINTER                   :: per_pot

      CALL timeset(routineN, handle)
      NULLIFY (grid, pw)
      dr1c = cgrid%pw_grid%dr(1)
      dr2c = cgrid%pw_grid%dr(2)
      dr3c = cgrid%pw_grid%dr(3)
      gbo = cgrid%pw_grid%bounds
      bo = cgrid%pw_grid%bounds_local
      grid2 => cgrid%array
      IF (par_scheme == do_par_atom) myind = 0
      Radius: DO IRadTyp = 1, SIZE(pgfs)
         per_pot => per_potentials(IRadTyp)%pot
         pw => per_pot%TabLR
         npts = pw%pw_grid%npts
         dr1 = pw%pw_grid%dr(1)
         dr2 = pw%pw_grid%dr(2)
         dr3 = pw%pw_grid%dr(3)
         grid => pw%array(:, :, :)
         !$OMP PARALLEL DO DEFAULT(NONE) &
         !$OMP SHARED(bo, gbo, grid, grid2, pw, npts, per_pot, mm_atom_index) &
         !$OMP SHARED(dr1, dr2, dr3, dr1c, dr2c, dr3c, par_scheme, mm_charges, mm_particles) &
         !$OMP SHARED(mm_cell, dOmmOqm, shells, para_env, IRadTyp, qmmm_spherical_cutoff) &
         !$OMP PRIVATE(Imm, LIndMM, IndMM, qt, sph_chrg_factor, ra, myind) &
         !$OMP PRIVATE(rt1, rt2, rt3, k, vec, ivec, xd1, xd2, xd3, ik1, ik2, ik3, ik4) &
         !$OMP PRIVATE(ij1, ij2, ij3, ij4, ii1, ii2, ii3, ii4, my_k, my_j, xs1, xs2, xs3) &
         !$OMP PRIVATE(p1, p2, p3, q1, q2, q3, r1, r2, r3, v1, v2, v3, v4, e1, e2, e3) &
         !$OMP PRIVATE(f1, f2, f3, g1, g2, g3, h1, h2, h3, s1, s2, s3, s4, a1, a2, a3) &
         !$OMP PRIVATE(b1, b2, b3, c1, c2, c3, d1, d2, d3, t1, t2, t3, t4, u1, u2, u3, val) &
         !$OMP PRIVATE(rv1, rv2, rv3, abc_X, abc_X_Y)
         Atoms: DO Imm = 1, SIZE(per_pot%mm_atom_index)
            IF (par_scheme == do_par_atom) THEN
               myind = Imm + (IRadTyp - 1)*SIZE(per_pot%mm_atom_index)
               IF (MOD(myind, para_env%num_pe) /= para_env%mepos) CYCLE Atoms
            END IF
            LIndMM = per_pot%mm_atom_index(Imm)
            IndMM = mm_atom_index(LIndMM)
            qt = mm_charges(LIndMM)
            IF (shells) THEN
               ra(:) = pbc(mm_particles(LIndMM)%r - dOmmOqm, mm_cell) + dOmmOqm
            ELSE
               ra(:) = pbc(mm_particles(IndMM)%r - dOmmOqm, mm_cell) + dOmmOqm
            END IF
            ! Possible Spherical Cutoff
            IF (qmmm_spherical_cutoff(1) > 0.0_dp) THEN
               CALL spherical_cutoff_factor(qmmm_spherical_cutoff, ra, sph_chrg_factor)
               qt = qt*sph_chrg_factor
            END IF
            IF (ABS(qt) <= EPSILON(0.0_dp)) CYCLE Atoms
            rt1 = ra(1)
            rt2 = ra(2)
            rt3 = ra(3)
            LoopOnGrid: DO k = bo(1, 3), bo(2, 3)
               my_k = k - gbo(1, 3)
               xs3 = REAL(my_k, dp)*dr3c
               my_j = bo(1, 2) - gbo(1, 2)
               xs2 = REAL(my_j, dp)*dr2c
               rv3 = rt3 - xs3
               vec(3) = rv3
               ivec(3) = FLOOR(vec(3)/pw%pw_grid%dr(3))
               xd3 = (vec(3)/dr3) - REAL(ivec(3), kind=dp)
               ik1 = MODULO(ivec(3) - 1, npts(3)) + 1
               ik2 = MODULO(ivec(3), npts(3)) + 1
               ik3 = MODULO(ivec(3) + 1, npts(3)) + 1
               ik4 = MODULO(ivec(3) + 2, npts(3)) + 1
               p1 = 3.0_dp + xd3
               p2 = p1*p1
               p3 = p2*p1
               q1 = 2.0_dp + xd3
               q2 = q1*q1
               q3 = q2*q1
               r1 = 1.0_dp + xd3
               r2 = r1*r1
               r3 = r2*r1
               u1 = xd3
               u2 = u1*u1
               u3 = u2*u1
               v1 = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*p1 + 12.0_dp*p2 - p3)
               v2 = -22.0_dp/3.0_dp + 10.0_dp*q1 - 4.0_dp*q2 + 0.5_dp*q3
               v3 = 2.0_dp/3.0_dp - 2.0_dp*r1 + 2.0_dp*r2 - 0.5_dp*r3
               v4 = 1.0_dp/6.0_dp*u3
               DO j = bo(1, 2), bo(2, 2)
                  xs1 = (bo(1, 1) - gbo(1, 1))*dr1c
                  rv2 = rt2 - xs2
                  vec(2) = rv2
                  ivec(2) = FLOOR(vec(2)/pw%pw_grid%dr(2))
                  xd2 = (vec(2)/dr2) - REAL(ivec(2), kind=dp)
                  ij1 = MODULO(ivec(2) - 1, npts(2)) + 1
                  ij2 = MODULO(ivec(2), npts(2)) + 1
                  ij3 = MODULO(ivec(2) + 1, npts(2)) + 1
                  ij4 = MODULO(ivec(2) + 2, npts(2)) + 1
                  e1 = 3.0_dp + xd2
                  e2 = e1*e1
                  e3 = e2*e1
                  f1 = 2.0_dp + xd2
                  f2 = f1*f1
                  f3 = f2*f1
                  g1 = 1.0_dp + xd2
                  g2 = g1*g1
                  g3 = g2*g1
                  h1 = xd2
                  h2 = h1*h1
                  h3 = h2*h1
                  s1 = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*e1 + 12.0_dp*e2 - e3)
                  s2 = -22.0_dp/3.0_dp + 10.0_dp*f1 - 4.0_dp*f2 + 0.5_dp*f3
                  s3 = 2.0_dp/3.0_dp - 2.0_dp*g1 + 2.0_dp*g2 - 0.5_dp*g3
                  s4 = 1.0_dp/6.0_dp*h3
                  DO i = bo(1, 1), bo(2, 1)
                     rv1 = rt1 - xs1
                     vec(1) = rv1
                     ivec(1) = FLOOR(vec(1)/pw%pw_grid%dr(1))
                     xd1 = (vec(1)/dr1) - REAL(ivec(1), kind=dp)
                     ii1 = MODULO(ivec(1) - 1, npts(1)) + 1
                     ii2 = MODULO(ivec(1), npts(1)) + 1
                     ii3 = MODULO(ivec(1) + 1, npts(1)) + 1
                     ii4 = MODULO(ivec(1) + 2, npts(1)) + 1
                     !
                     ! Spline Interpolation
                     !

                     a1 = 3.0_dp + xd1
                     a2 = a1*a1
                     a3 = a2*a1
                     b1 = 2.0_dp + xd1
                     b2 = b1*b1
                     b3 = b2*b1
                     c1 = 1.0_dp + xd1
                     c2 = c1*c1
                     c3 = c2*c1
                     d1 = xd1
                     d2 = d1*d1
                     d3 = d2*d1
                     t1 = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*a1 + 12.0_dp*a2 - a3)
                     t2 = -22.0_dp/3.0_dp + 10.0_dp*b1 - 4.0_dp*b2 + 0.5_dp*b3
                     t3 = 2.0_dp/3.0_dp - 2.0_dp*c1 + 2.0_dp*c2 - 0.5_dp*c3
                     t4 = 1.0_dp/6.0_dp*d3

                     abc_X(1, 1) = grid(ii1, ij1, ik1)*v1 + grid(ii1, ij1, ik2)*v2 + grid(ii1, ij1, ik3)*v3 + grid(ii1, ij1, ik4)*v4
                     abc_X(1, 2) = grid(ii1, ij2, ik1)*v1 + grid(ii1, ij2, ik2)*v2 + grid(ii1, ij2, ik3)*v3 + grid(ii1, ij2, ik4)*v4
                     abc_X(1, 3) = grid(ii1, ij3, ik1)*v1 + grid(ii1, ij3, ik2)*v2 + grid(ii1, ij3, ik3)*v3 + grid(ii1, ij3, ik4)*v4
                     abc_X(1, 4) = grid(ii1, ij4, ik1)*v1 + grid(ii1, ij4, ik2)*v2 + grid(ii1, ij4, ik3)*v3 + grid(ii1, ij4, ik4)*v4
                     abc_X(2, 1) = grid(ii2, ij1, ik1)*v1 + grid(ii2, ij1, ik2)*v2 + grid(ii2, ij1, ik3)*v3 + grid(ii2, ij1, ik4)*v4
                     abc_X(2, 2) = grid(ii2, ij2, ik1)*v1 + grid(ii2, ij2, ik2)*v2 + grid(ii2, ij2, ik3)*v3 + grid(ii2, ij2, ik4)*v4
                     abc_X(2, 3) = grid(ii2, ij3, ik1)*v1 + grid(ii2, ij3, ik2)*v2 + grid(ii2, ij3, ik3)*v3 + grid(ii2, ij3, ik4)*v4
                     abc_X(2, 4) = grid(ii2, ij4, ik1)*v1 + grid(ii2, ij4, ik2)*v2 + grid(ii2, ij4, ik3)*v3 + grid(ii2, ij4, ik4)*v4
                     abc_X(3, 1) = grid(ii3, ij1, ik1)*v1 + grid(ii3, ij1, ik2)*v2 + grid(ii3, ij1, ik3)*v3 + grid(ii3, ij1, ik4)*v4
                     abc_X(3, 2) = grid(ii3, ij2, ik1)*v1 + grid(ii3, ij2, ik2)*v2 + grid(ii3, ij2, ik3)*v3 + grid(ii3, ij2, ik4)*v4
                     abc_X(3, 3) = grid(ii3, ij3, ik1)*v1 + grid(ii3, ij3, ik2)*v2 + grid(ii3, ij3, ik3)*v3 + grid(ii3, ij3, ik4)*v4
                     abc_X(3, 4) = grid(ii3, ij4, ik1)*v1 + grid(ii3, ij4, ik2)*v2 + grid(ii3, ij4, ik3)*v3 + grid(ii3, ij4, ik4)*v4
                     abc_X(4, 1) = grid(ii4, ij1, ik1)*v1 + grid(ii4, ij1, ik2)*v2 + grid(ii4, ij1, ik3)*v3 + grid(ii4, ij1, ik4)*v4
                     abc_X(4, 2) = grid(ii4, ij2, ik1)*v1 + grid(ii4, ij2, ik2)*v2 + grid(ii4, ij2, ik3)*v3 + grid(ii4, ij2, ik4)*v4
                     abc_X(4, 3) = grid(ii4, ij3, ik1)*v1 + grid(ii4, ij3, ik2)*v2 + grid(ii4, ij3, ik3)*v3 + grid(ii4, ij3, ik4)*v4
                     abc_X(4, 4) = grid(ii4, ij4, ik1)*v1 + grid(ii4, ij4, ik2)*v2 + grid(ii4, ij4, ik3)*v3 + grid(ii4, ij4, ik4)*v4

                     abc_X_Y(1) = abc_X(1, 1)*t1 + abc_X(2, 1)*t2 + abc_X(3, 1)*t3 + abc_X(4, 1)*t4
                     abc_X_Y(2) = abc_X(1, 2)*t1 + abc_X(2, 2)*t2 + abc_X(3, 2)*t3 + abc_X(4, 2)*t4
                     abc_X_Y(3) = abc_X(1, 3)*t1 + abc_X(2, 3)*t2 + abc_X(3, 3)*t3 + abc_X(4, 3)*t4
                     abc_X_Y(4) = abc_X(1, 4)*t1 + abc_X(2, 4)*t2 + abc_X(3, 4)*t3 + abc_X(4, 4)*t4

                     val = abc_X_Y(1)*s1 + abc_X_Y(2)*s2 + abc_X_Y(3)*s3 + abc_X_Y(4)*s4
                     !$OMP ATOMIC
                     grid2(i, j, k) = grid2(i, j, k) - val*qt
                     !$OMP END ATOMIC
                     xs1 = xs1 + dr1c
                  END DO
                  xs2 = xs2 + dr2c
               END DO
            END DO LoopOnGrid
         END DO Atoms
         !$OMP END PARALLEL DO
      END DO Radius
      CALL timestop(handle)
   END SUBROUTINE qmmm_elec_with_gaussian_LG

! **************************************************************************************************
!> \brief Compute the QM/MM electrostatic Interaction collocating
!>      (1/R - Sum_NG Gaussians) on the coarser grid level.
!>      Long Range QM/MM Electrostatic Potential with Gaussian - Low Level
!> \param pgfs ...
!> \param grid ...
!> \param mm_charges ...
!> \param mm_atom_index ...
!> \param mm_particles ...
!> \param para_env ...
!> \param potentials ...
!> \param mm_cell ...
!> \param dOmmOqm ...
!> \param par_scheme ...
!> \param qmmm_spherical_cutoff ...
!> \param shells ...
!> \par History
!>      07.2004 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   SUBROUTINE qmmm_elec_with_gaussian_LR(pgfs, grid, mm_charges, mm_atom_index, &
                                         mm_particles, para_env, potentials, &
                                         mm_cell, dOmmOqm, par_scheme, qmmm_spherical_cutoff, shells)
      TYPE(qmmm_gaussian_p_type), DIMENSION(:), POINTER  :: pgfs
      TYPE(pw_r3d_rs_type), INTENT(IN)                   :: grid
      REAL(KIND=dp), DIMENSION(:), POINTER               :: mm_charges
      INTEGER, DIMENSION(:), POINTER                     :: mm_atom_index
      TYPE(particle_type), DIMENSION(:), POINTER         :: mm_particles
      TYPE(mp_para_env_type), POINTER                    :: para_env
      TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER       :: potentials
      TYPE(cell_type), POINTER                           :: mm_cell
      REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: dOmmOqm
      INTEGER, INTENT(IN)                                :: par_scheme
      REAL(KIND=dp), DIMENSION(2), INTENT(IN)            :: qmmm_spherical_cutoff
      LOGICAL                                            :: shells

      CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_elec_with_gaussian_LR'

      INTEGER                                            :: handle, i, Imm, IndMM, IRadTyp, ix, j, &
                                                            k, LIndMM, my_j, my_k, myind, n1, n2, &
                                                            n3
      INTEGER, DIMENSION(2, 3)                           :: bo, gbo
      REAL(KIND=dp)                                      :: dr1, dr2, dr3, dx, qt, r, r2, rt1, rt2, &
                                                            rt3, rv1, rv2, rv3, rx, rx2, rx3, &
                                                            sph_chrg_factor, Term, xs1, xs2, xs3
      REAL(KIND=dp), DIMENSION(3)                        :: ra
      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: pot0_2
      REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: grid2
      TYPE(qmmm_pot_type), POINTER                       :: pot

      CALL timeset(routineN, handle)
      n1 = grid%pw_grid%npts(1)
      n2 = grid%pw_grid%npts(2)
      n3 = grid%pw_grid%npts(3)
      dr1 = grid%pw_grid%dr(1)
      dr2 = grid%pw_grid%dr(2)
      dr3 = grid%pw_grid%dr(3)
      gbo = grid%pw_grid%bounds
      bo = grid%pw_grid%bounds_local
      grid2 => grid%array
      IF (par_scheme == do_par_atom) myind = 0
      Radius: DO IRadTyp = 1, SIZE(pgfs)
         pot => potentials(IRadTyp)%pot
         dx = Pot%dx
         pot0_2 => Pot%pot0_2
         !$OMP PARALLEL DO DEFAULT(NONE) &
         !$OMP SHARED(pot, par_scheme, para_env, mm_atom_index, mm_particles, dOmmOqm, mm_cell, qmmm_spherical_cutoff) &
         !$OMP SHARED(bo, gbo, dr1, dr2, dr3, grid2, shells, pot0_2, dx, mm_charges, IRadTyp) &
         !$OMP PRIVATE(myind, Imm, LIndMM, IndMM, ra, qt, sph_chrg_factor, rt1, rt2, rt3, my_k, my_j) &
         !$OMP PRIVATE(rv1, rv2, rv3, rx2, rx3, r, r2, rx, Term, xs1, xs2, xs3, i, j, k, ix)
         Atoms: DO Imm = 1, SIZE(pot%mm_atom_index)
            IF (par_scheme == do_par_atom) THEN
               myind = Imm + (IRadTyp - 1)*SIZE(pot%mm_atom_index)
               IF (MOD(myind, para_env%num_pe) /= para_env%mepos) CYCLE Atoms
            END IF
            LIndMM = pot%mm_atom_index(Imm)
            IndMM = mm_atom_index(LIndMM)
            ra(:) = pbc(mm_particles(IndMM)%r - dOmmOqm, mm_cell) + dOmmOqm
            qt = mm_charges(LIndMM)
            IF (shells) &
               ra(:) = pbc(mm_particles(LIndMM)%r - dOmmOqm, mm_cell) + dOmmOqm
            ! Possible Spherical Cutoff
            IF (qmmm_spherical_cutoff(1) > 0.0_dp) THEN
               CALL spherical_cutoff_factor(qmmm_spherical_cutoff, ra, sph_chrg_factor)
               qt = qt*sph_chrg_factor
            END IF
            IF (ABS(qt) <= EPSILON(0.0_dp)) CYCLE Atoms
            rt1 = ra(1)
            rt2 = ra(2)
            rt3 = ra(3)
            LoopOnGrid: DO k = bo(1, 3), bo(2, 3)
               my_k = k - gbo(1, 3)
               xs3 = REAL(my_k, dp)*dr3
               my_j = bo(1, 2) - gbo(1, 2)
               xs2 = REAL(my_j, dp)*dr2
               rv3 = rt3 - xs3
               DO j = bo(1, 2), bo(2, 2)
                  xs1 = (bo(1, 1) - gbo(1, 1))*dr1
                  rv2 = rt2 - xs2
                  DO i = bo(1, 1), bo(2, 1)
                     rv1 = rt1 - xs1
                     r2 = rv1*rv1 + rv2*rv2 + rv3*rv3
                     r = SQRT(r2)
                     ix = FLOOR(r/dx) + 1
                     rx = (r - REAL(ix - 1, dp)*dx)/dx
                     rx2 = rx*rx
                     rx3 = rx2*rx
                     Term = pot0_2(1, ix)*(1.0_dp - 3.0_dp*rx2 + 2.0_dp*rx3) &
                            + pot0_2(2, ix)*(rx - 2.0_dp*rx2 + rx3) &
                            + pot0_2(1, ix + 1)*(3.0_dp*rx2 - 2.0_dp*rx3) &
                            + pot0_2(2, ix + 1)*(-rx2 + rx3)
                     !$OMP ATOMIC
                     grid2(i, j, k) = grid2(i, j, k) - Term*qt
                     !$OMP END ATOMIC
                     xs1 = xs1 + dr1
                  END DO
                  xs2 = xs2 + dr2
               END DO
            END DO LoopOnGrid
         END DO Atoms
         !$OMP END PARALLEL DO
      END DO Radius
      CALL timestop(handle)
   END SUBROUTINE qmmm_elec_with_gaussian_LR

END MODULE qmmm_gpw_energy
